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A  method  of  formulating  and  automatically  integrating  the  equations 
of  motion  of  quite  general  constrained  dynamic  systems  is  presented. 

Design  sensitivity  analysis  is  also  carried  out  using  a  state  space  method 
that  has  been  used  extensively  in  structural  design  optimization.  Both 
dynamic  analysis  and  design  sensitivity  analysis  and  optimization  are 
shown  to  be  well-suited  to  application  of  efficient  sparse  matrix  computa¬ 
tional  methods.  Numerical  integration  is  carried  out  using  a  stiff 
numerical  integration  method  that  treats  mixed  systems  of  differential 
and  algebraic  equations.  A  computer  code  that  implements  the  method  for 
planar  systems  is  outlined  and  a  numerical  example  is  treated.  The  dynamic 
response  of  a  classical  slider-crank  is  analyzed  and  its  design  is 
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[.  Introduction  to  Computational  Dynamics 

Computational  methods  for  dynamic  analysis  of  electronic  circuits  [1] 
and  structures  [2]  have  become  well  developed  and  user-oriented,  a  situation 
that  is  not  enjoyed  in  the  field  of  mechanical  system  dynamics.  In  most 
areas  of  mechanical  design,  classical  Lagrangian  methods  of  dynamic  analysis 
are  still  used  almost  exclusively.  It  is  the  purpose  of  this  paper  to 
present  an  integrated  method  of  constrained  mechanical  system  dynamic 
analysis  and  design  sensitivity  analysis  that  is  user  oriented  and  capable 
of  routine  application  to  large  scale  systems. 

Development  of  general  purpose  computed  codes  for  dynamic  analysis 
of  mechanisms  and  mechanical  systems  are  initiated  in  the  mid-1960's  and 
has  made  significant  progress  for  certain  classes  of  mechanical  systems. 

A  comprehensive  survey  of  computer  codes  and  the  analytical  techniques  on 
which  they  are  based  has  been  presented  by  Paul  [3].  He  cites  numerous 
planar  (2-D)  mechanism  analysis  codes,  but  only  two  three-dimensional  (3-D) 
analysis  codes.  At  the  present  time  only  the  3-D  Integrated  Mechanism 
Program  (IMP)  [A]  and  the  2-D  Dynamic  Response  of  Articulated  Machinery 
(DRAM)  [5]  codes  are  well  developed  and  user-oriented.  However,  they  are 
both  based  on  closed  loop  linkages  or  mechanisms,  which  is  too  restrictive 
a  class  for  general  use  in  mechanical  system  dynamic  analysis  and  design. 

A  general  3-dimensional  Automatic  Dynamic  Analysis  of  Mechanical  Systems 
(ADAMS)  method  [6]  has  been  presented  that  is  fundamentally  better  suited 
for  dynamic  analysis  of  systems  that  are  not  in  closed  loop  configurations. 
However,  no  generally  applicable  computer  code  is  available  at  the 
present  time.  The  ADAMS  modelling  method  is  all  the  more  attractive 
since  it  is  well-suited  for  extension  to  design  sensitivity  analysis. 
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Before  presenting  the  theory  it  is  helpful  to  constrast  two  funda¬ 
mentally  different  approaches  to  modelling  linkages  and  mechanical  systems 
that  are  presently  in  use.  The  loop  closure  method  generates  equations 
that  require  closure  of  each  independent  loop  of  the  linkage.  The 
resulting  nonlinear  equations  are  then  differentiated  to  obtain  the 
smallest  possible  number  of  independent  equations  of  motion,  in  terms  of 
a  minimum  number  of  system  degrees  of  freedom.  Thus,  the  loop  closure 
method  generates  a  small  number  of  highly  nonlinear  equations  that  are 
solved  with  standard  numerical  integration  methods.  The  constrained 
system  modelling  method,  on  the  other  hand,  explicitly  treats  three 
degrees  of  freedom  for  each  element  of  the  system  (in  the  plane). 

Algebraic  equations  prescribing  constraints  between  the  various  bodies 
are  then  written  and  elementary  forms  of  equations  of  motion  for  each 
body  are  written  separately.  The  constraint  equations  prescribing 
assembly  of  the  mechanism  are  adjoined  to  the  equations  of  motion  through 
use  of  Lagrange  multipliers.  Thus,  one  treats  a  large  number  of  equations 
in  many  variables.  These  equations  may  be  solved  by  an  implicit  numerical 
integration  method  [1]  that  iteratively  solves  a  linear  matrix  equation. 
The  saving  grace  of  this  technique  is  that  the  matrix  that  arises  in  the 
iterative  solution  process  is  sparse.  That  is,  only  three  to  ten  percent 
of  the  elements  of  the  matrix  are  different  from  zero. 

It  is  the  purpose  of  this  paper  to  develop  the  theory  of  dynamic  and 
design  sensitivity  analysis  and  to  present  a  computer  code  called  Dynamic 
Analysis  and  Design  System  (DADS)  that  implements  the  constrained  system 
modelling  method  for  planar  systems.  The  same  basic  theory  is  applicable 
for  three-dimensional  systems,  but  due  to  analytical  complexity  it  will  be 
presented  in  a  separate  paper. 
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2.  Constrained  Equations  of  Planar  Motion 

For  planar  mechanical  systems,  constraints  between  elements  are  taken 
as  friction  free  (workless)  translational  and  rotational  joints.  Exten¬ 
sions  to  include  constraints  such  as  cams,  prescribed  functional  relations, 
and  intermittent  motion  are  possible  by  incorporating  provisions  for  non¬ 
standard  elements  that  are  supplied  by  the  user.  In  addition  to  standard 
constraints,  springs  and  dampers  connecting  any  pair  of  points  on  different 
bodies  of  the  system  are  included  in  the  model.  In  addition  to  these 
standard  force  elements,  allowance  is  made  for  external  forcing  functions. 

Generalized  Coordinates:  In  order  to  specify  the  configuration  or 
state  of  a  planar  mechanical  system,  it  is  first  necessary  to  define 
generalized  coordiates  that  specify  the  location  of  each  body.  As  shown 
in  Fig.  2.1,  let  the  x-y  coordinate  system  be  an  inertial  reference  frame. 
Define  a  body  fixed  coordinate  system  embedded  in  body  i,  with 

its  origin  at  the  center  of  mass.  One  can  now  locate  the  rigid  body  i  by 
specifying  the  global  coordinates  (x^.y^)  of  the  origin  of  the  body  fixed 
coordinates  and  the  angle  ip^  of  rotation  of  this  sytem  relative  to  the 
global  coordinates. 

Let  p  be  a  point  on  body  i,  as  shown  in  Fig.  2.1.  The  coordinates 
of  the  point  in  the  body  fixed  system  are  and  rj  The  same 

point  may  also  be  located  by  its  global  coordinates  x^  and  y^*.  The 
relation  between  these  coordinates  is 
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where  E(<j>^)  is  the  transformation  matrix 


E(4>±) 


cos  -  sin  4^' 


sin  cos 

l  i- 


(2.2) 


In  terms  of  these  generalized  coordinates,  one  can  write  the  kinetic 
energy  of  the  ic^  body  as 

1  *2*2  1  *2 

^i  *  2  “i(xi  +yi)+2Ji*i  (2,3) 

where  m^  is  the  mass  of  the  ic^  body  and  is  its  centroidal  polar 
moment  of  inertia. 

Further,  one  can  write  the  virtual  work  of  all  externally  applied 
forces  on  body  i  as 


6W  *0  5x.  +  Q  6y  +  Qx  <5 

e  i  y±  i  4^  *1 


(2.4) 


The  effect  of  all  forces  except  the  workless  constraint  forces  between 

bodies  is  included  on  <5W  of  Eq.  2.4. 

e  n 

Equations  of  Constraint:  Figure  2.2  depicts  two  adjacent  bodies  i 

and  j .  The  origins  of  their  body  fixed  coordinate  systems  are  located  by 

the  vectors  R^  and  R^  with  respect  to  the  inertial  frame  of  reference. 

Let  an  arbitrary  point  p^  on  body  i  be  located  by  r^  and  p^  on  body  j 

be  located  by  r...  These  points  are  in  turn  connected  by  a  vector  r  . 

ji  P 

One  can  write  a  vector  equation  beginning  at  the  origin  of  the  inertial 
reference  frame  and  closing  there,  to  obtain  the  vector  relationship 


Ri  +  rij  +  rP  -  rji  -  Rj  ‘  ° 


(2.5) 


*rr~T 


5 


The  constraint  equations  for  a  revolute  joint  are  now  obtained  by  re¬ 
quiring  that  and  coincide.  Setting  r^  =  0  and  writing  Eq.  2.5  in 
component  form,  one  has 


X.  +  5..  COS  <(>  -T)  sin  4>.  -  x.  -  cos  +  n..  sin  d>.  ■  0 

1  ij  iij  ijji  Jji  j 


y±  +  sin  cos  <j>i  -  y..  -  ^  sin  <f^  -  cos  $  *  0^ 


(2.6) 


For  a  translational  joint  shown  in  Fig.  2.3,  let  points  p  and  p„ 

lie  on  some  line  parallel  to  the  path  of  relative  motion  between  the  two 

bodies.  In  addition,  let  them  be  located  such  that  r. .  and  r..  are  per¬ 
il  Ji 

pendicular  to  this  line  and  of  nonzero  magnitude.  Successively  forming 
the  dot  product  of  Eq.  2.5  with  r^  and  r^  and  adding,  one  obtains  the 
scalar  equation 


rji  +  rji,<Ri-V 

(2.7) 


which  reduces  to 


cos  <j>i  -  Hy  sin  <(>i  +  cos  <(>..  -  sin  (xt  -  xj 
+  (5^  sin  <t>i  +  hij  cos  <^  +  5^  sin  4^+n^cos  >  (y±  -  y^) 


,  .  2  .  2  2  2  . 
+  Ei)  +  "u  '  -  V  ■ 0 


(2.8) 


A  second  scalar  equation  is  obtained  by  noting  the  r^  x  *  0.  Expan¬ 
sion  of  the  cross  product  yields  only  a  z  component,  which  must  be  zero. 


This  is 
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(^ij  cos  <p±  -  sin  sin  cos  $  ) 

-  (C^  cos  <jK  -  sin  sin  4>i  +  cos  )  -  0  (2.9) 


Spring-Dampers ;  Since  springs  and  dampers  generally  appear  in  pairs, 
they  are  incorporated  into  a  single  set  of  equations.  If  one  or  the  other 
is  absent,  its  effect  is  eliminated  by  setting  that  term  to  zero.  The 
equations  for  spring-damper  force  and  torque  are 


Fij  "  [kijUij  "  40ij)  +  cij  vij  +  F 


1^ 

ijJ  *1J 


(2.10) 


-  v  *  v«  * 


(2.11) 


where 


F  is  the  resultant  force  vector  [F  ,F  l1  in  the  spring- 


damper 


Xij  V 


R  is  the  vector  [£..  cos  a,  l. ,  sin  a]  between  points 
3ij  ij  ij 

and  of  a  spring-damper  connection  on  the  two 
bodies,  as  in  Fig.  2.4 

is  the  torque  acting  on  two  bodies  at  a  revolute  joint 

k  and  k  are  elastic  spring  coeffients 

2  rij 

c  and  c  are  damping  coefficients 

U  rij 

ln  and  q>n  are  the  undeformed  spring  length  and  angular 

ij  ij 

rotation  in  the  revolute 

i and  are  the  deformed  spring  length  and  angular 

rotation  in  the  revolute,  and  v  and  $  are  the  time 
derivatives  of  z  and  $ 
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F  and  T  are  constant  forces  and  torques  applied  along 

ij  ij 

the  spring  and  around  the  revolute  joint  between  two 
bodies 

From  Fig.  2.4  a  vector  expression  similar  to  Eq.  2.5  is  written  as 


R.  +  r  +  R  -r  -R.=0 
i  s,.  s . .  s . .  ; 

ij  IJ  Ji 


or  in  component  form 


Z,  .  cos  a 

ij 


Z . .  sin  a 
.  iJ 


x.  £  cos  4> .  -  n  sin  <p . 

IS  IS  1 

Ij  ij 


y.  |£  sin  <)>.  +  n  cos  <t. 
ij  Ls..  i  s . .  ij 

13  13 


cos  4  . 

-  n  sin  4. 

‘ji 

3 

Sji  3 

sin  4 , 

+  n  cos  4. 

'ji 

3 

sji  J 

(2. 12' 


where  a  is  the  angle  between  R  and  the  inertial  x  axis.  Equation  2.12 

Sij 

is  used  to  obtain  Z..  by  noting  that 
il 

2  2  i/2 

l..  =  [(&. .  cos  a)  +  (£..  sin  a)  ] 
ij  ij  ij 

=  [  (-x .  -  £  cos  4 .  +  n  sin  <p .  +  x.  +  £  cos  <p . 
i  s..  i  s..  i  i  s..  i 

ij  il  3i 


-  n  sin  <p  )  +  (-  y  -  £  sin  <j>  -  n  cos  <p 

Sji  3  ij  1  Sij  X 

2  1/2 

+  y.  +  £  sin  4.  +  n  cos  A.)  ] 

3  Sji  j  Sji  3 


(2.13) 


Substituting  the  left  side  of  Eq.  2.12  into  Eq.  2.10,  the  following  force 
expressions  are  obtained: 


F*  *  !F-mI  cos  a  (2.14) 

ij  J 


F 

y 


ij 


I  sin  ® 


(2.15) 


where  cos  a  and  sin  a  are  obtained  by  dividing  Eq. 
defining 


2.12  by  l... 

ij 


Finally, 


ij 


*ij 


(2.16) 


and  transfering  F 

ij  x 


,  and  v..  to  the  right  hand  sides  of  Eqs. 


ij  yij 

2.13  to  2.16,  one  obtains  equations  in  the  form  required  by  the  numerical 
integration  algorithm. 


System  Equations  of  Motion:  With  the  kinetic  energy  given  by  Eq. 

2.3  and  the  generalized  forces  associated  with  externally  applied  forces 

and  spring-dampers,  one  can  construct  the  equations  of  motion  for  the 

system.  To  allow  for  a  unified  development,  let  q^  denote  the  vector 
T 

^Xi,yi’^i^  oF  generalized  coordinates  of  body  i.  Denote  the  kinetic 
energy  of  body  i  as  KE^q*),  the  generalized  force  on  body  i  as  Qi(qi), 
and  the  equations  of  constraint  between  bodies  i  and  j  as  ‘^(q^’.q^)  *  0, 
i  <  j.  Here  the  constraint  numbering  convention  is  that  i  <  j  in  order 
to  preclude  inclusion  of  the  same  set  of  constraint  equations  twice. 

Presuming  that  the  constraints  are  workless,  one  has  the  variational 


equation  of  motion  [7] 
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3KE , 


l  Tt  (tt)  'o1-  ?  ^‘V-o 

3q 


(2.17) 


that  must  hold  for  all  time  and  for  all  virtual  displacements  dq  that  are 
consistent  with  the  constraints.  This  is 


,  i  ,  3«lj 


r  dq"1-  +  ~~  -  dqJ  =  0  ,  i  <  j 

3q  3  q^ 


(2.18) 


By  the  Farkas  Lemmas  of  optimization  theory  [9],  there  exist  multipliers 
Alj ,  i  <  j ,  such  that 


,  /3KE,  \  .  .T  .  , 

V#r  -v* +  v 


.T 


J 


3$  J  .  i  .  3*  -  .  i 
- r  <5q  + t  fiqJ 

3q  3qJ 


=  0  (2.19) 


for  all  dq*.  By  Eq.  2.3,  KEi  *  -j  q1  Miqi,  where 


M1  = 


m.  0  0 

0  m.  0 

l 

0  0  J. 


and  Eq.  2.19  may  be  written  in  the  form 


.T  .  .T 
q  M  -  Q1  +  l  x 


.  .T  -.ij  ,.T  a.ji 

ij  3$  +  £  ji  3$ 

i<j  V  j<i  X  Jq1 


dq1  =  0 


(2.20) 


Since  Eq.  2.20  must  hold  for  arbitrary  dq  ,  this  yields  the  equations 
of  motion 


mV  -  Q1  (q1)  +  l  ^—7  A1J  +  l  AJ1  -  0  ,  i-1 . r  (2.21) 

i<j  3q  j <i  3  q 


Jj  ^  V  a±£i  Ji 


These  differential  equations  and  the  algebraic  constraint  equations 


«ij(qi,qj)  -  o  ,  i  <  j 


(2.22) 


form  the  system  equations  of  motion.  Defining  the  generalized  velocity 
vector 


i 


u 


•i 

q 


one  can  replace  the  second  order  system  of  Eqs.  2.21  by  the  first  order 
system 


Mi.i 
M  u 


.1.  i. 

Q  (q  ) 


♦  l  ^  x1J 

i<j  3  q 


*  l  ^ 

j<i  3  q 


0 


(2.23) 


•i  i 
q  -  u 


(2.24) 


The  problem  of  dynamic  analysis  is  thus  reduced  to  solving  the 
system  of  differential  and  algebraic  equations  of  Eqs.  2.22  to  2.24,  with 
initial  conditions  given  that  are  consistent  with  the  constraints  of  Eqs. 
2.22.  The  solution  variables  are  q^  and  u^,  i«l,...,r  and  X*^  for  all 
i  <  j.  Each  of  these  variables  is  a  function  of  time. 
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3.  Integration  of  Constrained  Equations  of  Motion 

Numerical  Integration  of  Differential  Equations:  In  order  to  solve 
the  differential  equations  of  motion,  numerical  integration  theory  must 
be  used  to  obtain  a  set  of  approximations  that  are  suitable  for  digital 
computation.  Integration  of  the  equations  of  motion  must  be  accomplished 
as  asimultaneous  solution  of  the  algebraic  and  differential  equations  of 
Eqs.  2.22  to  2.24.  Standard  numerical  methods  are  however  designed  to 
solve  only  systems  of  differential  equations  of  the  form 

y  =  f(y,t)  (3.1) 

where  y  is  an  m-vector  of  unknowns  and  f  is  an  m-vector  of  known  functions. 
A  modified  approach  is  taken  here  that  allows  for  the  simultaneous  solution 
of  mixed  algebraic  and  differential  equations  of  the  form 

g(y,y,t)  =  0  (3.2) 

where  some  components  of  y  may  not  appear  in  any  of  the  equations,  or  they 
may  appear  nonlinearly.  Before  introducing  the  method  to  be  used,  it  is 
instructive  to  review  the  approach  applied  in  solving  Eq.  3.1. 

The  basic  method  of  constructing  approximate  solutions  [1]  is  to 
place  a  time  grid  t^,  i=l,...  on  the  interval  [0,t],  where  h^  =  t^+^  -  t^ 
is  the  grid  spacing.  One  then  approximates  the  solution  y(t)  of  Eq.  3.1 
on  the  time  grid  as  yi%y(t^). 

The  basic  approximating  equation  for  stiff  differential  equations  [1], 
i.e.  differential  equations  with  widely  split  eigenvalues,  is  the  Gear 
algorithm 
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i+1  .  a  *1-1-1  p 

y  -  i^Sq  y  -  l  « 
j-1  J 


i-j+l 


(3.3) 


where  Che  constants  8q  and  a ,  j*l,...,  k,  called  Gear's  coefficients,  are 
determined  so  that  Eq.  3.3  is  exact  for  any  polynomial  solution  of  Eq.  3.1 
of  degree  up  to  k.  Gear's  coefficients  [1]  also  have  the  property  that 
the  algorithm  tends  to  be  stable,  even  for  stiff  differential  equations. 

The  multistep  formula  that  is  used  to  solve  Eq.  3.2  is  derived  from 
Eq.  3.3.  One  progresses  from  t^  to  t^+^  by  solving  Eq.  3.3,  together  with 


,  i+1  .i+1  .  . 

g(y  ,  y  ,  ti+1)  -  0 


Linearization  of  Eq.  3.4  generates  the  Newton  formula 


te(m)  .  (m)  3g(m)  .(■)  (m) 

37“  4y  +  -*r—  Ayv  '  -  -  gw 


(3.4) 


(3.5) 


,  »  (m)  (m+1)  (m) 

where  Ay  *  y  -  y  and  m  represents  the  iteration  number.  The 

time  step  index  i  has  been  suppressed  here  for  notational  simplicity.  Sub¬ 
stitution  of  Eq.  3.5  into  Eq.  3.3,  noting  that  the  summation  term  of  Eq. 

3.3  remains  constant  at  each  iteration,  yields  the  corrector  formulas 


heQ  3y  3y 


..(m)  1  .  (m) 

4y  -  TT-  Ay  v 

heo 


(3.6) 


(3.7) 


If  Eq.  3.2  is  of  the  form  Py  +  f(y,t)  -  0,  then  Eq.  3.6  is  of  the 
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i  ’  ♦  *£]»«  -  -  .« 


(3.8) 


The  iterative  corrector  procedure  is  continued  at  each  time  step  until 
all  of  the  Newton  differences  Ay^  are  below  a  specified  tolerance  level. 
At  each  iteration  y  and  y  are  updated  as 


r<«H-l)  =  (m)  +  (a) 


.(m+1)  ,  •<»)  +  A-(m)  ,  -(a)  +  Ay(a) 

h60 


It  is  interesting  to  note  that  every  element  of  y  is  obtained  at  each 
time  step,  even  though  many  elements  of  y  may  not  appear  in  any  of  the 
equations. 


Sparse  Matrix  Algebra:  The  sytem  of  nonlinear  algebraic  and  differ¬ 
ential  equations  defined  in  the  previous  paragraphs  is  very  loosely 
coupled.  For  two  reasons,  no  attempt  is  made  to  eliminate  variables  to 
obtain  a  smaller  system  of  equations.  First,  the  Jacobian  matrix  that  is 
formed  by  linearization  of  the  loosely  coupled  equations  and  used  in 
iterative  solution  by  Newton’s  method  is  sparse  and  can  be  very  efficiently 
stored  and  decomposed.  Second,  the  repetitive  nature  of  the  loosely 
coupled  equations  results  in  compact  and  efficient  computer  routines  for 
evaluating  the  nonlinear  equations  and  nonzero  matrix  entries.  Recently 
developed  sparse  matrix  algorithms  [8]  make  both  of  these  operations 
practical  and  desirable.  It  has  been  shown  [8]  that  it  is  usually  more 
efficient  to  solve  large  systems  of  sparse  equations,  rather  than  smaller 
systems  with  greater  percentages  of  nonzero  entries. 


Consideration  of  matrix  sparsity  is  important  to  the  speed  of  computa¬ 
tion  in  problems  of  dynamic  system  analysis  [1].  When  less  than  30%  of  the 
matrix  entries  are  nonzero,  it  is  inefficient  to  store  the  matrix  as  a  full 
two-dimensional  array.  Instead,  only  the  nonzero  entries  are  stored  in 
compacted  form.  The  method  most  commonly  used  for  compacting  the  data  is 
to  store  the  row  and  column  indices  of  each  nonzero-valued  entry  in  the 
matrix  in  two  vectors  I  and  J  and  the  value  in  a  third  vector  A.  This  is 
called  "i-J"  ordering. 

Sparse  matrix  algorithms  are  most  efficient  when  the  nonzero  matrix 
entries  are  stored  in  an  organized  manner.  This  usually  implies  that  they 
be  evaluated  row  by  row  or  column  by  column.  The  previously  mentioned 
repetitive  matrix  evaluation  scheme  is  defeated  by  this  requirement,  since 
it  usually  results  in  the  evaluation  of  small  submatrices  located  at 
various  positions  throughout  the  matrix.  To  overcome  this  difficulty,  a 
permutation  vector  is  generated  from  the  row  and  column  vectors  describing 
the  nonzero-valued  positions.  As  each  matrix  entry  is  generated,  it  is 
directed  to  a  specific  location  in  the  "A"  vector  by  a  permutation  index. 

At  completion  of  the  matrix  evaluation,  all  entries  are  stored  exactly 
as  if  they  had  been  evaluated  in  column  order. 

A  sparse  matrix  decomposition  algorithm  is  then  applied  to  the  column 
ordered  matrix  and  an  LU  factorization  [1]  is  accomplished.  Full  pivoting 
is  not  achieved,  but  the  algorithm  chooses,  from  among  the  largest  of  the 
acceptable  pivot  elements,  the  pivot  that  results  in  a  minimum  number  of 
fills  in  the  resulting  L  and  U  matrices.  This  is  important  for  efficient 
execution  of  the  forward  and  backward  substitution  phases,  since  an  in¬ 
creased  number  of  fills  destroys  the  original  matrix  sparsity  and  results 
in  excessive  computer  time. 
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4.  Elements  of  Design  Sensitivity  Analysis  and  Optimization 

As  shown  In  this  section  and  developed  in  more  detail  in  companion 

papers  [9],  the  constrained  dynamic  system  model  of  Section  2  and  the 

implicit  numerical  integration  formulation  of  Section  3  are  ideally  suited 

for  efficient  design  sensitivity  analysis.  If  the  designer  identifies  a 

set  of  parameters  b  =  [b^,...,b  ]  that  he  uses  to  specify  the  system 

design,  a  natural  question  that  arises  is:  How  does  a  design  variation 
T 

6b  *  [6b^,...,6b  ]  effect  the  dynamic  response  of  the  system?  An  even 
more  important  question  is:  What  is  the  effect  on  dynamic  response  of  a 
design  variation  6b  that  is  required  to  be  consistent  with  certain  per¬ 
formance  requirements?  A  method  of  design  sensitivity  analysis  developed 
for  related  classes  of  problems  [10]  is  employed  here  to  answer  these 
questions. 

In  order  to  address  the  design  sensitivity  analysis  problem,  a  compact 

T 

notation  is  required.  Let  the  vectors  z(t)  =  [z,  (t) , . . . ,  z  (t)]  contain 

1  m 

all  generalized  displacement  and  velocity  coordinates; 


£(t)  =  [ £^(t) , ^(t) » • • • 1  contain  spring  lengths,  velocities  and  forces; 

i  -t  T  j 

and  A(t)  *  [A  J  (t),  i  <  j]  contain  Lagrange  multipliers  that  are 


associated  with  constraints  (which  determine  reaction  forces  associated 
with  constraints). 

With  this  notation,  the  equations  of  motion  are  written  in  compact 
form  as 


P(b)z  +  f (t,z,A,i,b)  *  0 
z(0)  -  v(b) 


(4.1) 


The  matrix  P  is  made  up  of  6  x  6  matrices  that  are  associated  with 
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derivative  terms  in  Eqs.  2.23  and  2.24  for  each  of  the  bodies;  i.e., 
P^b)  -  diag  (m^(b) ,  ra^b),  J^b),  1,  1,  1).  Similarly,  the  vector 
is  made  up  of  the  remaining  terms  in  Eqs.  2.23  and  2.24;  i.e.. 


(Q„  -  I 


3<t 


ij 


i<j  3xi 


-  I 


3$' 


ji 


j<i  3xi 


xji)  . 


<Q„  -  l 


m13  . 


1<J  3yl 


j<l  3yi 


%  -  I 

*i  i<j 


34> 


ij 


3$. 


t  ij 


3*Ji  , ji. 

-r—  X  )  ,  u  ,  v  ,  w 
j<i  3*i  i  i  i 


-  I 


The  spring  damper  equations  of  Eqs.  2.13  to  2.16  in  matrix  form  are 


Ii  +  9 (z, £,b)  -  0 

Ci  ■  vi(b)-  j  -  ° 


•  1 >  *  -  * » s-l  | 


(4.2) 


where  s  is  the  number  of  spring  damper  pairs  and  the  matrix  I  is  diagonal 
and  contains  all  zeros,  except  for  ones  on  the  diagonal  corresponding  to 
coefficients  of  in  Eq.  2.16.  The  vector  0  simply  contains  all  other 
terms  of  Eqs.  2.13  to  2.16. 

Finally,  the  set  of  kinematic  constraints  is  written  in  vector  form  as 


$(z,b)  -  0  (4.3) 

The  dependence  of  Eqs.  4.1  to  4.3  on  the  design  variable  b  represents  the 
fact  that  b  may  include  such  parameters  as  spring  and  joint  locations, 

i 

spring  and  damper  coefficients,  masses,  and  moments  of  inertia. 

Optimal  Design  Problem  Formulation:  For  the  purpose  of  this  intro¬ 


ductory  discussion  of  design  sensitivity  analysis,  the  cost  function  to 


be  minimized  and  performance  constraints  to  be  satisfied  are  taken  in 


integral  form.  The  general  design  problem  is  then  to  minimize 


80(b) 


■T 

■  n 


1*Q  (  t  ,  Z  ,  A  ,  Z  ,  b )  d  t 


(4.4) 


subject  to  the  constraints 


♦i 


g^b)  + 


T 

L  (t,z, A,  it,b)dt 

J0 


(  m  0*  i=1» • • • »r' 

|  <  0,  i-r'+l, . . . ,r 


(4.5) 


This  form  of  integral  constraint  can  be  made  to  include  constraints 
that  must  hold  for  all  time;  i.e.. 


n(t,z  (t),A(t),£(t),b)  <  0  ,  0  <  t  <  t  (4.6) 

The  inequality  of  Eq.  4.6  is  equivalent  to  the  integral  constraint 

•T 

[  n(t,z(t)  ,  A(t)  ,  Z(t)  ,b)  +  |n(t,z(t),X(t),£.(t),b)  |]dt  =  0 

J0 

which  is  of  the  form  of  Eq.  4.5.  This  technique  has  been  used  for  a  wide 
variety  of  constrained  dynamic  systems  [11]  with  good  success.  It  may  be 
noted  that  constraints  of  the  form  of  Eq.  4.6  may  represent  excursion 
limits,  bounds  on  stress,  path  generation  error  bounds,  and  a  multitude 
of  other  forms  of  meaningful  design  constraints. 

Design  Sensitivity  Analysis:  In  order  to  determine  the  effect  of  a 
design  variation  6b,  one  must  first  note  that  the  solution  of  Eqs.  4.1 
to  4.3  will  be  changed  by  perturbations  6z,  5i,,  and  51.  Presuming  the 
design  problem  is  well  posed,  small  perturbations  6b  in  b  will  lead  to 
small  perturbations  in  the  solution  variables. 
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Consider  a  typical  functional  of  Eq.  4.4  or  4.5  where  the  subscript 
has  been  deleted  for  notational  simplicity.  To  first  order,  one  has 

s*.j£ib+f  [f  +  f  «  +  ft  «  +  f  «-]  <*■» 


The  objective  is  to  eliminate  explicit  dependence  of  Eq.  4.7  on  &z,  dA, 

and  61  and  to  write  explicitly  in  terms  of  db.  To  do  this,  an  effective 

—  * 

technique  is  to  introduce  adjoint  variables  p,  u  ,  and  p  ,  multiplying 
Eqs.  4.1  to  4.3  respectively,  and  Integrate  to  obtain  the  identities 


ft 


{pT[P(b)z  +  f  (t,z, A,2,b)) }dt  *  0 


[  {p,T[H  +  0(z,A,b)]}dt  -  0 
J0 

[  {y*T$(z,b) }dt  *  0 
JO 


(4.8) 


Then  integrate  Eqs.  4.8  by  parts  and  require  that  u,  u',  and  p*  satisfy 
the  adjoint  equations 


T  - 

-P  u  + 


3fT  - 


a. 


+  4* 


p*  + 


39 


3L 


0 


3_£  -  3LT 
3z  U  ~  3A 


0 


(4.9) 


and  terminal  conditions 
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y(x)  =  0,u'  (t)  =  0,  y*(i)  =  0 
The  funcational  5iJ<  may  then  be  rewritten  as 

<$i(i  =  5b 


(4. 10) 


(4.11) 


where 


t  +  ;T(0)p(b)  J*  +  Y  u;+4.(0)  ^±-1 


9b 


fT 

*9L 

-T  9(P(b)z) 

-T  3f 

„*T  li  - 

,T  30 

0 

l3b ' 

U  9b 

”  y  3b  ~ 

U  3b 

M  3b 

dt  (4.12) 


One  may  now  integrate  the  differential  equations  4.9  from  r  to  0, 
using  the  terminal  conditions  of  Eq.  4.10.  If  one  uses  the  implicit 
numerical  algorithm  of  Eq.  3.8  for  the  above  set  of  equations,  the  Newton 
iteration  formula  is 


ZdLpT  +*±T) 

T 

ii 

T1 

3# 

Au 

\h8Q  3z  / 

3z 

dz 

3fT 

(zL.  tt  +  iiT\ 

r\ 

Ay’ 

3 1 

lh6o  3*  / 

U 

dZ 

0 

0 

Ay* 

(4.13) 


where  g  represents  the  left  hand  side  of  Eq.  4.9.  For  ready  reference,  it 
is  desirable  to  repeat  Eq.  3.8  applied  to  Eqs.  4.1  to  4.3  here  as 


where  g  represents  the  left-hand  side  of  Eqs.  4.1  to  4.3. 


Note  that  the  coefficient  matrix  in  Eq.  4.13  is  exactly  the  transpose 
of  the  matrix  of  Eq.  4.14  with  the  exception  of  the  (-1/hBg)  term.  But 

since  one  integrates  from  t  to  0,  this  implies  negative  time  steps  h,  hence 

3  f  -  3Q 

one  need  not  modify  the  expressions  (l/hg„  P  +  — )  and  (1/hB^  I  +  -— )  of 

0  3z  0  3  2. 

Eq.  4.14  for  use  in  Eq.  4.13.  Thus,  one  can  use  the  same  sparse  matrix 
factorization  code  generated  for  dynamic  analysis  to  solve  Eq.  4.13  during 


adjoint  calculations.  The  numerical  efficiencies  are  obvious.  This  entire 
subject  is  addressed  in  detail  in  Ref.  9. 

Method  of  Generalized  Steepest  Descent:  The  linearized  optimal 
design  problem  is  to  find  a  vector  6b  that  minimizes 


6tj>g  -A  5b 


subject  to  design  constraints 


6ip  ■  A  6b 


i  -  -  'I'g.  6-1,. ...r' 

(<  -  1 1>  ,  Be(r'+1, . . .  ,r) ,  <l>  >  - 


and  the  quadratic  step-size  constraint 


6bTW6b  <  C2 


(4.15) 


(4.16) 


(4.17) 
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where  <p  is  a  vector  of  e-active  constraints,  A  contains  sensitivity  coef¬ 
ficients  associated  with  e-active  constraints,  W  is  a  positive  definite 
weighting  matrix,  and  £  and  e  are  small  positive  numbers. 

Employing  the  Kuhn-Tucker  necessary  conditions  of  nonlinear  program¬ 
ming  [11],  one  derives  the  necessary  design  change  6b  as 


6b 


+  6b 


(4.18) 


where 


6b1  =  W"1(A°  +  A  y1) 


2  -1-2 
5b  =  -W  A  y 


„  1  -1  .0 
M  Y  =  -A  W  A 
W 


M  y  ==  -  <  ji 

w 


'T  -1  - 
M,  ,  =  A  W  A 


(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 


and  - —  >  0  is  to 

2yo 

design  derivative 
errors. 


be  chosen  as  a  step-size.  The  vector  6b1-  is  a  constrained 
2 

and  6b  is  a  vector  design  change  that  corrects  constraint 


\ 
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5.  Dynamic  Analysis  and  Design  System  (DADS) 

The  DADS  computer  program  implements  the  dynamic  analysis,  sensitivity 
analysis,  and  optimal  design  methods  presented  in  Sections  2  to  4. 

Figures  5.1  and  5.2  are  diagrams  showing  the  subprograms  that  are  used 
for  dynamic  and  sensitivity  analysis.  Figure  5.3  gives  the  overall  program 
flow  diagram  that  incorporates  these  subprograms. 

The  dynamic  analysis  phase  (DYNANL)  of  the  program  generates  sparse 
matrix  code  for  pivoting  and  LU  factorization  and  solves  the  system  of 
differential  -  algebraic  equations  for  the  state  variables  during  a 
specified  time  interval.  It  employs  sparse  matrix  codes  and  the  numerical 
integration  algorithm  of  Section  3.  In  addition,  the  Jacobian  matrix  and 
state  variables  at  each  time  step  are  stored  on  a  direct  access  disk  file 
for  subsequent  use  in  sensitivity  analysis. 

The  sensitivity  analysis  and  optimization  phase  t S5NANL >  .\‘  the  pro¬ 
gram  carries  out  the  calculations  of  Section  4.  Thu  program  solves  the 
system  of  adjoint  equations  using  the  previously  stored  data  and  the  same 
numerical  integration  routine  as  above.  The  program  further  computes  the 
sensitivity  coefficients  and  necessary  design  improvements.  This  process 
is  repeated  until  an  optimum  design  is  achieved. 

The  Dynamic  Analysis  Phase:  The  DYNANL  phase  of  DADS  establishes  the 
sparse  matrix  code  description  of  the  mechanical  system  and  numerically 
solves  the  differential  and  algebraic  equations  for  the  state  variables. 

As  shown  in  Fig.  5.1,  this  involves  two  major  steps:  (i)  generation  of 
an  initial  sparse  code  description  (including  pivoting  and  LU  factoriza¬ 
tion  code)  and  (ii)  repetitive  solution  of  the  Newton  iteration  equations 
for  the  state  variables  during  the  time  interval  of  interest. 


\ 
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In  Che  first  step  of  DYNANL,  estimates  of  the  initial  configuration 
of  the  system  are  provided  by  INDATA  and  used  by  VARSET  to  initialize  a 
state  variable  vector  for  subsequent  use  by  the  numerical  integration 
routine  DIFSUB.  A  compact  numbering  system  identifying  bodies,  joints, 
and  spring-damper  elements  is  used  to  input  data  through  INDATA  and  pro¬ 
vides  the  necessary  description  of  the  mechanical  system  configuration. 
This  information  is  used  to  construct  the  Newton  corrector  equations,  Eqs. 
4.14  and  adjoint  equations,  Eqs.  4.13,  for  use  in  the  sensitivity  analsyis 
phase.  The  components  of  these  equations  are  obtained  from  the  equations 
of  motion,  Eqs.  2.23  and  2.24,  spring-damper  equations,  Eqs.  2.11  and  2.13 
to  2.16,  and  constraint  equations,  Eqs.  2.7  to  2.9.  The  program  then  uses 
subroutine  S03000  to  generate  initial  vectors  of  row  and  column  indices 
that  locate  nonzero  entries  in  the  Jacobian  matrix.  Similarly,  user 
supplied  row-column  positions  of  nonstandard  elements  are  provided  by 
incorporating  the  necessary  code  in  USET.  A  symbolic  description  of  the 
resultant  matrix  is  printed  by  DEBUGG  for  reference  purposes  and  a  column 
ordering  permutation  vector  is  generated  in  S08000  (Section  3). 

Subroutine  S01000  evaluates  coefficients  and  the  right  hand  side  of 
Eq.  3.8,  or  equivalently  Eq.  4.14.  Its  functions  are  to  (i)  evaluate 
force  and  displacement  functions  of  time  that  are  provided  by  the  user 
(FOREXT) ,  (ii)  transfer  the  state  variables  (Section  2)  from  a  single 
vector  used  by  the  numerical  integration  routine  to  the  standard  variables 
and  user-supplied  variables  (US0LV1) ,  (iii)  evaluate  the  Jacobian  matrix 
for  the  standard  equations  and  user-supplied  equations  (US0LV2)  using 
updated  variables  from  step  (ii)  and  the  previously  generated  (S08000) 
column  ordering  permutation  vector,  and  (iv)  evaluate  the  standard 
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equations  and  user-supplied  equations  (US0LV3).  Finally  a  sparse  LU 
factored  description  of  the  matrix  is  generated  in  INVERT. 

The  second  step  of  DYNANL  is  to  numerically  integrate  the  system 
of  equations  during  the  time  interval  of  interest.  This  is  accomplished 
by  the  numerical  integration  subroutine  DIFSUB,  which  repeatedly  calls 
S01000  to  update  the  Jacobian  matrix  and  system  of  equations  as  it 
executes  the  iterative  corrector  formula  of  Eq.  4.14. 

The  Sensitivity  Analysis  Phase:  In  the  SENANL  phase,  DYNANL  is 
again  called  for  adjoint  calculation  by  AHLPSI,  where  the  design  sensi¬ 
tivity  coefficients  are  calculated.  For  this  case,  DYNANL  reads  the 
previously  generated  data  from  disk  file  and  executes  its  two  step 
procedure.  In  the  first  step,  sparse  matrix  codes  for  pivoting  and  LU 
factorization  are  generated  for  the  transpose  of  the  Jacobian  matrix, 
using  INVERT.  In  the  second  step,  the  system  of  adjoint  equations  is 
numerically  integrated  by  DIFSUB.  The  program  repeatedly  reads  data 
from  the  disk,  calls  S01NEW  for  reevaluation  of  the  right  hand  side  of 
the  adjoint  system  of  equations  and  iteratively  solves  for  the  adjoint 
variables  at  each  time  step.  The  adjoint  variables  are  then  stored  on 
disk  for  later  use  in  calculation  of  the  sensitivity  coefficients. 

Description  of  the  DADS  Program:  The  Dynamic  Analysis  and  Design 
System  Incorporates  the  two  phases  of  Sections  5. 2  and  5.3  into  a  design 
program  that  is  capable  of  handling  a  variety  of  planar  design  problems. 

A  brief  description  of  how  these  phases  are  coupled  as  a  design  tool  is 
given  below.  For  more  detailed  descriptions,  reference  9  is  recommended. 
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Initial  estimates  and  bounds  on  design  paramters,  the  number  of 
constraints  and  other  parameters  related  to  the  design  problem  (Section  4) 
are  read  by  the  main  program.  Then  system  parameters  such  as  masses, 
moments  of  inertia,  location  of  centers  of  mass,  applied  constant  forces, 
joint  types,  and  spring-damper  parameters  are  read  through  subroutine 
INDATA.  Subroutine  RELATE  relates  the  variables  or  other  system  parameters 
of  the  dynamic  analysis  phase  DYNANL  (Section  2)  to  the  updated  design 
parameters.  For  example,  if  the  location  of  a  revolute  joint  on  a  link 
is  changed,  it  may  be  necessary  to  change  its  n.-ir.:-.  and  moment  of  inertia. 
This  step  is  required  following  each  design  iteration.  In  the  subsequent 
dynamic  analysis  phase  the  state  equations  are  solved  by  DYNANL  and  the 
Jacobian  matrix,  state  variables,  time  step,  current  time,  and  order  of 
the  numerical  integration  algorithm  are  stored  on  a  direct  access  disk. 

This  information  is  first  used  to  evaluate  the  integrands  of  the  cost  and 
constraint  functionals  of  Eqs.  4.4  and  4.5  and  the  integral  form  of  Eq. 

4.6  through  subroutine  HALS. 

The  inequality  functional  constraints  of  Eqs.  4.5  and  4.5  are  then 

tested  for  violation  and  the  corresponding  adjoint  equations  of  Eqs.  4.9 

and  4.10  are  solved  (TEQFUN,  INFUNC,  ADJONT).  The  sensitivity  coefficient 

matrix  of  Eq.  4.12  is  calculated  for  each  e-active  functional  constraint 

(AHLPSI).  Subroutine  DPARMC  tests  each  of  the  design  parameter  constraint 

and  calculates  the  corresponding  sensitivity  matrices.  The  sensitivity 

vector  of  Eq.  4.15  is  evaluated  by  AHLPSI.  Subroutine  NEWB  computes 

1  2 

the  matrix  M  of  Eq.  4.23,  solves  Eq.  4.21  and  4.22  for  y  and  y  ,  and 
computes  the  necessary  design  changes  6b  given  by  Eqs.  4.18  to  4.20. 
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At  this  stage  the  convergence  criteria  are  tested.  If  they  are 
satisfied,  the  new  design  is  taken  as  the  optimum.  Otherwise,  the  process 
is  repeated  with  the  new  design  parameters  as  the  initial  estimate. 

Brief  descriptions  are  given  for  other  important  subprograms  appearing 
in  Fig.  5.2: 

AHSM:  The  function  AHSM  determines  the  maxima  of  the  integrands  of  the 
transformed  equality  constraints  obtained  from  Eq.  4.6.  These  maxima 
are  used  in  TEQFUN. 

FUNPSE:  The  function  FUNPSE  evaluates  the  cost  and  constraint  functionals 
of  Eqs.  4.4  and  4.5. 

PBFN:  The  subroutine  PBFN  calculates  terms  of  the  matrix  P(b)  appearing 
in  Eq.  4.1. 

ADPDB:  The  subroutine  ADPDB  calculates  derivatives  of  the  elements  of 
P(b)  with  respect  to  the  design  parameters  b. 

DGDBZ:  The  subroutine  DGDBZ  calculates  the  derivatives  of  the  non  integral 
parts  of  the  cost  and  constraint  functionals  of  Eqs.  4.4  and  4.5. 

DNUDB:  The  subroutine  DNUDB  is  used  only  when  initial  conditions  of  the 
dynamic  analysis  are  dependent  on  design  parameters.  It  evaluates  the 
derivatives  of  the  initial  values  of  the  solution  variables  with  respect  to 
design  parameters. 

3L.  3L 

ADLDZ:  The  subroutine  ADLDZ  calculates  the  derivatives  —  ,  —  ,  and 

o  Z  o  A 


6.  Applications  and  Numerical  Results 

The  DADS  program  has  been  developed  to  treat  analysis  and  design  of 
quite  general  planar  dynamic  systems.  To  test  the  program  a  relatively  simple 
slider  crack  mechanism  is  analyzed,  design  sensitivity  analysis  is  carried 
out  and  the  design  is  optimized. 

The  radial  slider-crank  mechanism  is  a  four-bar  linkage  of  rigid 
bodies  that  move  in  a  plane.  Fig.  6.1  shows  the  approximate  initial 
position  of  such  a  mechanism  with  one  spring-damper  pair.  Link  1  is 
ground,  link  2  is  the  crank  shaft,  link  3  is  the  connecting  rod  or 
coupler,  and  link  4  is  the  piston  or  slider.  A  spring-damper  pair  is 
attached  between  link  4  and  ground,  as  shown  in  Fig.  6.1.  Revolute  joints 
connect  bodies  1  and  2,  2  and  3,  3  and  4.  A  translational  joint  connects 
bodies  4  and  1.  Gravitational  forces  are  ignored  in  the  present  simulation. 

A  symbolic  listing  of  the  nonzero  positions  of  the  Jacobian  matrix 
for  this  example  is  given  in  Fig.  6.2.  The  only  significance  to  be 
attached  to  the  digits  and  letters  is  that  there  is  a  nonzero  entry  in 
each  of  the  noted  positions.  All  other  entries  are  zeros,  so  only  9.5% 
of  the  matrix  elements  are  nonzero  and  are  accounted  for  by  sparse  matrix 
methods. 

Formulation  of  the  Optimal  Design  Problem:  By  virtue  of  its  move¬ 
ment,  a  radial  slider-crank  mechanism  exerts  a  force  on  ground  through 
the  crank  bearing  and  the  wrist-pin  guide.  It  is  desirable  to  keep  these 
"shaking  forces"  within  bounds.  It  is  also  desired  that  an  upper  bound 


be  placed  on  the  angular  velocity  of  the  crank  at  the  final  instant  t  of 
the  tine- interval  [0,t]  under  consideration.  The  cost  function  is  chosen 


co  be  twice  the  maximum  energy  that  is  stored  in  the  spring  during  the 
interval  of  motion.  The  design  parameters  shown  in  Fig.  6.1,  are  as 
follows: 
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-  spring  constant  k^  of  the  spring 
b^  *  height  of  the  points  of  attachment  of  the  spring 
b^  *  half  of  the  length  of  the  coupler 
With  the  notation  of  Sections  2  and  4,  the  optimal  design  problem  is 
stated  as  follows:  Minimize 


*0  5  ma*  V*1  -  lioy 

0  <  t  <  T 


(6.1) 


subject  to  the  equations  of  motion,  Eqs.  2.23  and  2.24,  the  equations  of 
constraint,  Eqs.  2.6  to  2.9,  the  spring-damper  relations,  Eqs.  2.13  to 
2.16,  the  initial  conditions  of  the  form  of  the  second  equation  of  Eqs. 
4.2,  the  functional  constraints 


|A*2|  <  A*2(«ax), 


0  <  t  <  T 


(6.2) 


4>2(t)  <  *2  (max) 


(6.3) 


and  the  design  parameter  constraints 


bi  -  bi  -  bi  * 


1  -  1,2,3 


(6.4) 


Here  is  the  underformed  spring  length,  is  the  current  spring  length, 

and  is  the  angular  velocity  of  the  crankshaft.  From  Eqs.  2.6  and  2.13 

1  12 

3  $  "12  12 

the  expression  — —  A^  ■  A^  is  interpreted  as  the  y-component  of  the 
^  12 

generalized  force  at  this  joint.  Therefore,  A,,  is  the  y-component  of 
reaction  force  on  the  crankshaft  at  the  joint  between  the  crankshaft 
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and  ground.  For  the  sake  of  simplicity,  only  the  constraint  on  the 

vertical  component  of  shaking  force  at  the  crank  bearing  is  considered 
here. 

After  introduction  of  an  artificial  design  parameter  [11],  and  the 
integral  functional  forms  in  the  conventional  way  [11],  the  problem  can  be 
reformulated  as  follows:  Minimize 

*0  “  b4  (6-5> 

subject  to  Eqs.  2.23  to  2.24,  2.6  to  2.9,  2.13  to  2.16,  and  the  second 
equations  of  Eqs.  4.1  and  4.2;  the  constraints 


<l> 


1 


T  12  12 

<  -  A  -  X  (max)  >  dt  =  0 

0 


(6.6) 


<  X 


12 


-  A^Onax) 


>  dt 


(6.7) 


ft 


<  -  110>  -  b4  >  dt  -  0 


(6.8) 


I 

(6.9) 


i 

of  Eq.  6.4.  Here,  the  symbol  > 

!i 

4 


i* 


4»4  =  “  4>2^max^  -  0 

and  the  design  parameter  constraints 
<n(t)>  has  the  following  meaning: 
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<n(t)>  - 


(n(t))  ,  for  n(t)  >  0 

0  ,  for  n(t)  <  0 


(6.10) 


This  formulation  is  a  special  case  of  the  general  formulation  of  the 
optimal  design  problem  stated  in  Section  4.  After  specification  of  initial 
numerical  data,  the  DAOS  program  described  in  Section  5  can  be  implemented 
to  obtain  its  solution. 

Design  Sensitivity  Analysis:  In  optimal  design,  calculation  of 
design  derivatives  constitutes  a  principal  part  of  the  work  to  be  done. 

When  design  sensitivity  coefficients  are  obtained,  the  gradient  projection 
iterative  optimization  algorithm  of  Section  4  can  be  applied.  In  the 
present  problem,  subroutines  RELATE,  DNUDB,  DGDBZ,  ADLDZ,  and  DLDFB  are 
the  major  user-supplied  subprograms  that  are  required  for  design  sensi¬ 
tivity  analysis.  Among  them,  DGDBZ  and  ADLDZ  are  often  simple,  because 
they  depend  solely  on  the  forms  of  the  functional  constraints.  Sub¬ 
routines  DNUDB,  DGDBZ,  and  DLDFB  are  used  in  AHLPSI  to  calculate  sensi¬ 
tivity  coefficients. 

Numerical  Results:  Initial  estimates  of  the  parameters  for  the 
slider  crank  mechanism  under  consideration  are  given  in  Table  6.1.  The 
table  associates  the  computer  code  names  with  the  corresponding  parameters 
in  Eqs.  2.4,  2.6  to  2.9,  2.13  to  2.16,  and  2.23.  Two  simulation  time 
intervals  are  considered;  [0,1]  and  [0,2]  seconds.  During  transient 
analysis  a  constant  counter-clockwise  torque  of  100  in-lbf  is  applied 
to  link  2. 

To  study  the  numerical  results  of  design  sensitivity  analysis,  the 
sensitivity  coefficients  ,  Se(l,2,3,4),  and  l  ®  were  first  calculated 
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for  the  2  second  interval,  with  initial  estimates  of  the  design  parameters 

given  in  Table  6.2a.  Lower  and  upper  bounds  on  the  first  three  design 

L  T  U 

parameters  were  selected  as  b  =  [0.8,  0.2,  5.0]  and  b  =  [1.5,  0.8,  16.0]. 

12 

Also  A ^  (max)  =  10.0  and  ^^(raax)  =  0.3. 

The  cost  and  constraint  functional  were  evaluated  for  0.1%  and  1.0% 

perturbations  of  the  first,  second,  and  fourth  design  parameters  and  0.01% 

and  0.1%  changes  in  the  third  design  parameter.  The  corresponding  changes 

in  the  functionals  given  in  Table  6.2a  are  presented  in  compact  form.  In 

these  computations  the  integrands  of  \p and  have  been  normalized 
12  12  2 

by  dividing  by  (max),  (max),  and  b^,  respectively.  Sensitivity 
coefficients  are  given  in  Table  6.2b. 

The  predicted  changes  6if/  =  £  6b  in  the  cost  and  constraint  functions 

p 

due  to  the  foregoing  design  perturbations  were  calculated  and  are  given 

Ipo 

in  Table  6.2a.  It  is  observed  from  the  table  that  terms  of  the  form  £  6b 
match  satisfactorily  with  the  actual  changes  in  8  =0,3,4.  The  slight 

discrepancies  are  attributed  to  various  approximations  and  linearizations 
at  many  steps  of  both  transient  and  sensitivity  analysis  and  the  coarse 
time  grid  used  for  numerical  integration. 

In  carrying  out  the  optimal  design  algorithm  for  the  two  second 
interval,  a  design  reduction  ratio  of  3  percent  was  used  to  compute  the 
step  size  in  the  first  iteration.  The  initial  estimate  b°  and  design 
variable  bounds  are  given  in  Table  6.3. 

At  the  starting  design  ||6b"'"||  =  0.1248  and  |  1 6b2 1  |  =  0.2307.  Con¬ 
straints  3  and  4  are  violated  and  the  lower  bound  constraint  on  design 
parameter  b^  is  tight.  After  6  iterations,  as  the  algorithm  approached 


\a 

T;'  .  ■> W" : 
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1  _  C  O 

Che  optimum  design,  ||$bx||  reduced  to  0.8404  x  10  ,  ||6b  |j  reduced  to 

0.2669  x  10  \  and  design  variable  constraints  on  b^,  b ^  and  b^  became 

tight.  The  optimum  design  is  given  in  Table  6.3. 

Next,  the  design  problem  is  considered  for  the  interval  [0,1], 

with  design  and  initial  data  given  in  Table  6.4.  For  this  case  a  design 

reduction  ratio  .of  10Z  was  used  in  computing  the  step  size  in  the  first 

iteration.  At  the  starting  design  (|5b^||  =■  0.0,  ||6b^||  ■  0.1106  x  10 

and  constraint  3  is  violated.  In  the  second  iteration  | j 6b1 j  |  “  0.7348, 

2  _2 

| | 6b  j|  -  3.814  x  10  ,  and  again  constraint  3  is  violated.  After  18 

iterations,  as  the  algorithm  approaches  the  optimum  design,  ||6b^||> 

5.567  x  10  ,  | | 6b  ||  *  3.018  x  10  ,  and  design  parameter  constraints  on 

b^  and  b^  are  tight.  The  optimum  design  is  given  in  Table  6.4. 
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TABLE  6.1  Initial  Estimates  of  the  Parameters  for  the  Slider- 
Crank  Mechanism  (Units  in  Inch,  Pound-Force,  Second 
System) . 


LINK  DESCRIPTION 


=M(1) 

=1. 

m2  =M(2)  =4. 

^  =M(3)  =1. 

m,  =M(4) 

4 

=3. 

Jx  =JIN(1) 

=1. 

J2  =JIN(2)  =10. 

J3  =J1N(3)  =4. 

J.  =JIN(4) 

4 

=2. 

xx  =X(1) 

=0. 

x2  =X(2)  =1. 

x3  =X(3)  =8.667 

x4  =X(4) 

=17.32 

y1  =Y(D 

=0. 

y2  =Y(2)  =0. 

y3  =Y(3)  =5.25 

y4  =Y(4) 

=  .5 

^  =PHI(1) 

=0. 

<f2  =PHI  (2)  =1.5708 

<t>3  =PHI  (3)  =-.5236 

d>4  =PHI  (4) 

=0. 

QXl  =FX(1) 

=0. 

QX2  =FX(2)  =0. 

QX3  =FX(3)  =0. 

Qx^=FX(4) 

=0. 

Qyi  -FY(1) 

=0. 

Qy2  =FY(2)  =0. 

Qy  3  =FY(3)  =0. 

Qy4=FY(4) 

=0. 

%  -TQU> 

=0. 

Q  4>2  =TQ(2)  =100. 

=TQ(3)  =0. 

Q<b  =TQ(4) 

4 

=0. 

JOINT  DESCRIPTION 

i  -IB(1,1) 

=1 

i  =IB(1 , 2)=2 

i  =IB(1,3)=3 

i  =IB(1 ,4) 

=4 

c12  =X1(1) 

=0. 

C23  =X1(2)  =9.0 

=X1  (3)  =10. 

CA1“X1(4) 

=0. 

n12  -Y  1(1) 

=0. 

n23  =Y1(2)  =0. 

n34  =Y1 (3)  =0. 

n41=Yl(4) 

=0.5 

j  = IB (2 , 1) 

=2 

j  =IB(2,2)=3 

j  =IB(2,3)=4 

j  =IB(2,4) 

=1 

C2J_  =X2(1) 

=-l. 

C32  =X2(2)  =-10. 

C43  =X2(3)  =0. 

?i4=X2(4) 

=0. 

n21  =Y2  (1) 

=0. 

n32  =Y2(2)  =0. 

n43  =Y2 (3)  =0. 

V=Y2(4) 

1.0 

SPRING-DAMPER 

DESCRIPTION 

i  =1BSD(1 , 

D=1 

C-  =XF1 (1)  =35. 

s14 

nSl4=YFl(i)  =0.5 

j  =IBSD(2 , 

1)=4 

CS41=XF2(1) 

=2. 

%i=YF2a}  =o- 

k14  =SK(1)  =1. 

c14=dc(D 

=1. 

*14  -SDL(l) 

=15.68 

1  =SDL0(1)=15.68 

°14 

Fn  =SDF(1)  =0. 

U14 

*  The  parameters  FX,  FY,  and  TQ  represent  components  of  force  and  moment 
that  are  reduced  to  the  centers  of  mass  on  each  link  and  contribute  to 
QXi>  Qy^ >  and  Qij>^  in  2.4. 


a  ■if1*--  *— 


TABLE  6.2a  (continued) 


TABLE  6.2b  Sensitivity  Coefficients  l  and  l 


Components 


DYNANL 


SENANL 


Figure  5.3  Flow  Diagram  of  the  DADS  Computer  Program. 


Approximate  Initial  Configuration  of  the  Slider-Crank  Mechanism 
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Figure  6.2  Symbolic  Listing  of  the  Nonzero  Entries  in  the  Jacobian 
Matrix  for  the  Example  Slider-Crank  Mechanism 
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